library(RSocrata)
library(tidyverse)
library(lubridate)
library(tidygeocoder)
library(sf)
library(tmap)
library(spatstat)
source("Funciones_propias_ppp.R", encoding = "UTF-8")

Contextualización

La accidentalidad es uno de los problemas más comunes en el planeta provocando perdidas que pueden ir desde materiales hasta la vida misma, por lo que es de particular interés tratar de comprender como se comporta este fenómeno para evitar ser victimas de él.

La problemática en cuestión tiene la complicación de ser algo de naturaleza aleatoria, sin embargo esto no quiere decir que no exista alguna forma de controlarla pues la experiencia ha mostrado que las diferentes ciudades tienen lugares con mayor concentración de accidentes que otras ubicaciones dentro de la misma ciudad.

Debido a que se está trabajando la ocurrencia de accidentes en una ciudad de interés, la distribución Poisson y los procesos puntuales Poisson resultan ser una herramienta adecuada para tratar de explicar el fenómeno.

Obtención de la base de datos

Una vez presentado el fenómeno que se trabajará, se selecciona la ciudad de Barranquilla para realizar el estudio.

Los datos se extrajeron de la página web GOV.CO la cual tiene diferentes bases da datos de Colombia abiertas al público. ( Accidentalidad en Barranquilla). La base de datos contiene una gran cantidad de variables, sin embargo solo se escogen aquellas de interés para el patrón puntual. La estructura de los datos es presentada a continuación:

accidentes <- read.socrata("https://www.datos.gov.co/resource/yb9r-2dsi.json") %>%
    mutate(fecha = ymd(fecha_accidente)) %>%
    select(12, 6:8) %>%
    filter(year(fecha) == 2021)
knitr::kable(head(accidentes, 10),
             col.names = c("Fecha", "Gravedad del accidentes",
                           "Clase de accidentes", "Sitio accidente"))
Fecha Gravedad del accidentes Clase de accidentes Sitio accidente
2021-01-01 Solo daños Choque CL 31 CR 22
2021-01-01 Con heridos Choque CL 80B CR 42
2021-01-01 Con heridos Choque AV DEL RIO ALETA DEL TIBURON
2021-01-02 Con muertos Choque CL 74 21B 23
2021-01-02 Con heridos Choque CL 72 CR 50
2021-01-02 Solo daños Choque CR 18 CL 36B
2021-01-02 Con muertos Choque CL 23 CR 15
2021-01-02 Con heridos Choque VIA 40 CL 85
2021-01-02 Solo daños Choque CL 72 FRENTE AL 57 68
2021-01-03 Solo daños Choque CL 45 CR 4

Geocodificación de las direcciones

Se hace necesario realizarle unos ajustes a la base de datos considerada puesto que está no incluye las ubicaciones exactas de los accidentes ya sea en longitud - latitud o proyección UTM, para dicho propósito se usa la función geo() del paquete tidygeocoder para convertir las direcciones de los accidentes en coordenadas de longitud - latitud.

# NO CORRER ESTE CHUNK
direcciones <- paste(accidentes$sitio_exacto_accidente, ", Barranquilla, Colombia", sep = "")
localizaciones <- geo(address = direcciones, method = "arcgis")
write.csv(x = localizaciones, file = "accidentes.csv", row.names = F)

Luego de realizar obtener las coordenadas en longitud - latitud, se usan múltiples funciones del paquete sf para obtener las respectivas ubicaciones en formato UTM obteniéndose lo siguiente:

bqlla <- st_read('bqlla.geojson')%>%
  st_transform(crs = 3857)

coordenadas <- read.csv("accidentes.csv")

# Se toman 1as primeras 4700 filas para no tener que regenerar el csv
# de localizaciones
accidentes_sf <- cbind(accidentes, coordenadas[8927:13626, 2:3]) %>%
  st_as_sf(coords = c('long', 'lat')) %>%
  st_set_crs(value = 4326) %>%
  st_transform(crs = 3857) %>%
  st_intersection(bqlla) %>% 
  filter(gravedad_accidente == "Con heridos") %>% 
  select(-fecha, -gravedad_accidente)
knitr::kable(head(accidentes_sf[, c(1:2, 8)], 10))
clase_accidente sitio_exacto_accidente Shape_Area geometry
Choque CL 80B CR 42 146434812 POINT (-8329215 1230965)
Choque AV DEL RIO ALETA DEL TIBURON 146434812 POINT (-8323586 1231025)
Choque CL 72 CR 50 146434812 POINT (-8327065 1231853)
Choque VIA 40 CL 85 146434812 POINT (-8327354 1235362)
Choque CL 64 CR 66 146434812 POINT (-8325561 1232219)
Choque CR 38 CL 108 146434812 POINT (-8330597 1230153)
Choque CL 82 CR 59B 146434812 POINT (-8327880 1233398)
Choque AV CIUDAD CARIBE CL 125 146434812 POINT (-8330059 1230143)
Choque CL 30 CR 6 146434812 POINT (-8323371 1228407)
Choque CL 110 CR 9G 146434812 POINT (-8330826 1227455)

Gráfico de las localizaciones

tmap_mode('view')

tm_shape(bqlla)+
  tm_polygons(alpha = 0.3, border.alpha = 0.7)+
  tm_shape(shp = accidentes_sf)+
  tm_dots(size = 0.01)

Pruebas de homogeneidad en la intensidad

Uno de los parámetros críticos al momento de modelar un patrón puntual es la intensidad la cual puede ser constante (homogénea) o variable (inhomogénea), por lo tanto es importante iniciar el análisis con una prueba de homogeneidad de la intensidad.

# Definiendo el patron puntual de los datos
datos_ppp <- ppp(x = st_coordinates(accidentes_sf)[, 1],
             y = st_coordinates(accidentes_sf)[, 2],
             window = as.owin(W = bqlla))

Argumento gráfico

Para iniciar, se divide el mapa de Barranquilla en 25 cuadrantes, si el patrón fuera homogéneo se esperaría que el número de accidentes en cada uno de las subsecciones de la ciudad contenga aproximadamente la misma cantidad de accidentes.

ncuadrantes(datos_ppp)
##    x y
## 17 7 2
qc_datos <- quadratcount(datos_ppp, nx = 7, ny = 2)
plot(qc_datos, main = "Accidentes en 2021")

Se aprecia que existe discrepancia entre el número de accidentes por sectores. Claramente al este de la ciudad el número de accidentes es mayor respecto al oeste lo cual es señal de que el patrón puntual en cuestión es de naturaleza inhomogénea.

Prueba \(\chi^2\)

Adicional a los conteos del número de accidentes en cada uno de los sectores definidos previamente, se usa la prueba \(\chi^2\) para contrastar:

\[ \begin{cases} H_0: \text{La intensidad es constante} \\ H_1: \text{La intensidad no es constante} \end{cases} \]

A continuación se muestra el resultado obtenido

quadrat.test(qc_datos)
## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  
## X2 = 638.74, df = 12, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 
## Quadrats: 13 tiles (irregular windows)

Puesto que se tiene un valor - p demasiado pequeño (orden de \(10^{-16}\)) se rechaza la hipótesis de homogeneidad de intensidad, es decir, la intensidad puede ser modelada mediante una relación funcional \(\lambda(x,y)\) y se hace necesario estimarla con algún método ya sea paramétrico o no paramétrico como se vera más adelante.

Estimación de propiedades de primer orden

Previamente se verifico que el número de accidentes en la ciudad tiene intensidad inhomogénea, en esta sección el propósito es estimar dicha función ya sea con metodología paramétrica y no paramética

No paramétrica

A continuación se presentan mapas de probabilidad de accidentalidad al realizar una estimación no paramétrica de la función de intensidad basado en funciones kernel los cuales son de la forma:

\[ \hat{\lambda}(x) = \frac{1}{h^2} \sum_{i = 1}^{n} \frac{\kappa\left(\frac{||x-x_i||}{h}\right)}{q(||x||)} \]

El ancho de banda se seleccionó utilizando la función bw.scott, la cual está implementada para los problemas que involucra la estimación espacial de la intensidad de puntos observados en redes lineales, como es el caso de los accidentes viales. Además, se analizó el error estándar de las estimaciones de densidad obtenido con esta función comparado con los resultados obtenidos de otras posibles funciones disponibles en el paquete spatstat y se encontró que el menor error estándar promedio fue el hallado con la función seleccionada.

graph_ppp(datos_ppp, round(bw.scott(datos_ppp, isotropic = T), 2),
          "Año 2021")

Modelos log-lineales

Otra posible solución para estimar la función de intensidad es mediante modelos log - lineales los cuales son de la forma

\[ log \ \lambda (u) = \theta \cdot S(u) \]

donde \(S(u)\) es una función vectorial de valor real evaluada en alguna ubicación \(u = (x, y)\).

Para ajustar estos modelos en R, se usan la función ppm() del paquete spatstat la cual recibe como primer argumento un objeto de la clase ppp y como segundo una formula.

El ajuste de dichos es modelos es sensible a la escala de medición de las ubicaciones; debido a que las ubicaciones se encuentran en metros (los cuales son números muy grandes) es necesario cambiar la escala de medición por una más “pequeña” para evitar problemas de singularidad matricial por lo que se convierten los metros en kilómetros.

# Rescalando los objetos ppp para ajustar los modelos
datos_ppp_km <- rescale(datos_ppp, 1000)
unitname(datos_ppp_km) <- c("km", "kilometers")

Luego de reescalar las ubicaciones se ajustan algunos modelos log-lineales.

# Ajuste de modelos
# Tendencia lineal en x e y
log_linx <- ppm(datos_ppp_km, ~ x)
log_liny <- ppm(datos_ppp_km, ~ y)
log_lin_xy <- ppm(datos_ppp_km, ~ x + y)

## Tendencia cuadratica
log_lin_x_cuad <- ppm(datos_ppp_km, ~ x + I(x^2))
log_lin_y_cuad <- ppm(datos_ppp_km, ~ y + I(y^2))
log_lin_xy_cuad <- ppm(datos_ppp_km, ~ x + I(x^2) + y + I(y^2))

# Tendencia cuadratica en x pero no en y
log_lin_xcuad_y <- ppm(datos_ppp_km, ~ x + I(x^2) + y)
log_lin_inter <- ppm(datos_ppp_km, ~ x*y)
log_lin_cuad_inter <- ppm(datos_ppp_km,  ~ x*y + I(x^2) + I(y^2))

Lineal en \(x\)

\[ \log\ \lambda(x,y) = \beta_0 + \beta_1x \]

par(mfrow = c(1, 2)) 
plot(predict(log_linx), main = "Número de accidentes")
plot(datos_ppp_km, add = T, cex = 0.5)
diagnose.ppm(log_linx, which = "smooth", 
             main = "Residuales suavizados")

Lineal en \(y\)

\[ \log\ \lambda(x,y) = \beta_0 + \beta_2 y \]

par(mfrow = c(1, 2)) 
plot(predict(log_liny), main = "Número de accidentes")
plot(datos_ppp_km, add = T, cex = 0.5)
diagnose.ppm(log_liny, which = "smooth", 
             main = "Residuales suavizados")

Lineal en \(x\) e \(y\)

\[ \log\ \lambda(x,y) = \beta_0 + \beta_1 x + \beta_2 y \]

par(mfrow = c(1, 2)) 
plot(predict(log_lin_xy), main = "Número de accidentes")
plot(datos_ppp_km, add = T, cex = 0.5)
diagnose.ppm(log_lin_xy, which = "smooth", 
             main = "Residuales suavizados")

Cuadrático en \(x\)

\[ \log\ \lambda(x,y) = \beta_0 + \beta_1 x + \beta_{11} x^2 \]

par(mfrow = c(1, 2)) 
plot(predict(log_lin_x_cuad), main = "Número de accidentes")
plot(datos_ppp_km, add = T, cex = 0.5)
diagnose.ppm(log_lin_x_cuad, which = "smooth", 
             main = "Residuales suavizados")

Cuadrático en \(y\)

\[ \log\ \lambda(x,y) = \beta_0 + \beta_2 y + \beta_{22} y^2 \]

par(mfrow = c(1, 2)) 
plot(predict(log_lin_y_cuad), main = "Número de accidentes")
plot(datos_ppp_km, add = T, cex = 0.5)
diagnose.ppm(log_lin_y_cuad, which = "smooth", 
             main = "Residuales suavizados")

Cuadrático en \(x\) e \(y\)

\[ \log\ \lambda(x,y) = \beta_0 + \beta_1 x + \beta_2 y + \beta_{11} x^2 + \beta_{22} y^2 \]

par(mfrow = c(1, 2)) 
plot(predict(log_lin_xy_cuad), main = "Número de accidentes")
plot(datos_ppp_km, add = T, cex = 0.5)
diagnose.ppm(log_lin_xy_cuad, which = "smooth", 
             main = "Residuales suavizados")

Cuadrático en \(x\) lineal en \(y\)

\[ \log\ \lambda(x,y) = \beta_0 + \beta_1 x + \beta_{11} x^2 + \beta_2 y \]

par(mfrow = c(1, 2)) 
plot(predict(log_lin_xcuad_y), main = "Número de accidentes")
plot(datos_ppp_km, add = T, cex = 0.5)
diagnose.ppm(log_lin_xcuad_y, which = "smooth", 
             main = "Residuales suavizados")

Interacción entre \(x\) e \(y\)

\[ \log\ \lambda(x,y) = \beta_0 + \beta_1 x + \beta_2 y + \beta_{12} xy \]

par(mfrow = c(1, 2)) 
plot(predict(log_lin_inter), main = "Número de accidentes")
plot(datos_ppp_km, add = T, cex = 0.5)
diagnose.ppm(log_lin_inter, which = "smooth", 
             main = "Residuales suavizados")

Cuadrático con interacción en \(x\) e \(y\)

\[ \log\ \lambda(x,y) = \beta_0 + \beta_1 x + \beta_2 y + \beta_{12} xy+ \beta_{11} x^2 + \beta_{22} y^2 \]

par(mfrow = c(1, 2)) 
plot(predict(log_lin_cuad_inter), main = "Número de accidentes")
plot(datos_ppp_km, add = T, cex = 0.5)
diagnose.ppm(log_lin_cuad_inter, which = "smooth", 
             main = "Residuales suavizados")

Observaciones

Propiedades de segundo orden

Luego de realizar ajustes paramétricos y no paramétricos para la propiedad de primer orden (intensidad), ahora interesa estimar las propiedades de segundo orden las cuales son utiles para determinar la naturaleza del patrón puntual de interés: completamente aleatorio, inhibitorio o agregado.

# patron puntual no parametrico usado previamente
ppp_np <- density.ppp(datos_ppp, sigma = bw.scott,
                      isotropic = T)

Función K

# Implementacion manual de envolventes
# K inhomogena no parametrica
# Kmat_np <- matrix(0, nrow = 513, ncol = 99)
# radios <- seq(0, 3450, length.out = 513)
# 
# for(i in 1:99){
#   temp <- rpoispp(ppp_np)
#   temp <- Kinhom(temp, r = radios, lambda = ppp_np)
#   Kmat_np[, i] <- temp$bord.modif
# }
# Knp <- Kinhom(datos_ppp, lambda = ppp_np, r = radios)
# upperK <- apply(Kmat_np, 1, max)
# lowerK <- apply(Kmat_np, 1, min)
# 
# Knp_env_data <- data.frame(upper = upperK, 
#                            lower = lowerK,
#                            observed = Knp$bord.modif,
#                            theo = Knp$theo, 
#                            r = radios)
# saveRDS(Knp_env_data, "funciones_segundo_orden/Knp_env_data.Rds")

Knp_env_data <- readRDS("funciones_segundo_orden/Knp_env_data.Rds")

# Gráficando las curvas
# k_aux <- data.frame(r = rep(Knp_env_data$r, 2),
#                     values = c(Knp_env_data$theo,
#                                Knp_env_data$observed),
#                     kind = rep(c("Teórica", "Estimada"),
#                                each = length(Knp_env_data$r)))
# 
# Kgg <- ggplot(k_aux, aes(r, y = values, color = kind)) +
#   geom_ribbon(aes(ymin = rep(Knp_env_data$lower, 2),
#                   ymax = rep(Knp_env_data$upper, 2)),
#               color = "white", fill = "grey70") +
#   geom_path(aes(linetype = kind)) +
#   labs(x = "r [metros]",
#        y = "K", 
#        title = "Envolventes para la función K",
#        color = "" ,
#        linetype = "") +
#   theme_minimal() +
#   theme(plot.title = element_text(hjust = 0.5)) +
#   scale_color_manual(values = c("black", "red"))
# 
# Kggplotly <- plotly::ggplotly(Kgg)
# saveRDS(Kggplotly, "Kggplotly.Rds")

Kggplotly <- readRDS("Kggplotly.Rds")
Kggplotly

Función g

# Implementacion manual de envolventes
# g inhomogena no parametrica
# gnp <- pcfinhom(datos_ppp, lambda = ppp_np,
#                  r = radios)
# 
# gmat <- matrix(0, nrow = 513, ncol = 99)
# 
# for(i in 1:99){
#   temp <- rpoispp(ppp_np)
#   temp <- pcfinhom(temp, r = radios, lambda = ppp_np)
#   gmat[, i] <- temp$trans
# }
# 
# upperg <- apply(gmat, 1, max)
# lowerg <- apply(gmat, 1, min)
# 
# gnp_env_data <- data.frame(r = radios, 
#                            observed = gnp$trans,
#                            lower = lowerg, 
#                            upper = upperg, 
#                            theo = 1)
# saveRDS(gnp_env_data, "funciones_segundo_orden/gnp_env_data.Rds")

gnp_env_data <- readRDS("funciones_segundo_orden/gnp_env_data.Rds")

# Gráficando las curvas
# g_aux <- data.frame(r = rep(gnp_env_data$r, 2),
#                     values = c(gnp_env_data$theo,
#                                gnp_env_data$observed),
#                     kind = rep(c("Teórica", "Estimada"),
#                                each = length(gnp_env_data$r)))
# 
# g_gg <- ggplot(g_aux, aes(r, y = values, color = kind)) +
#   geom_ribbon(aes(ymin = rep(gnp_env_data$lower, 2),
#                   ymax = rep(gnp_env_data$upper, 2)),
#               color = "white", fill = "grey70") +
#   geom_path(aes(linetype = kind)) +
#   labs(x = "r [metros]",
#        y = "g", 
#        title = "Envolventes para la función g",
#        color = "" ,
#        linetype = "") +
#   ylim(c(0, 4.4)) +
#   theme_minimal() +
#   theme(plot.title = element_text(hjust = 0.5)) +
#   scale_color_manual(values = c("black", "red"))
# 
# g_plotly <- plotly::ggplotly(g_gg)
# saveRDS(g_plotly, "g_plotly.Rds")

g_plotly <- readRDS("g_plotly.Rds")
g_plotly